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Abstract 

We study the stability of a Bose condensate of atomic 7 Li in a (harmonic 
oscillator) magnetic trap at non-zero temperatures. In analogy to the stability 
criterion for a neutron star, we conjecture that the gas becomes unstable if 
the free energy as a function of the central density of the cloud has a local 
extremum which conserves the number of particles. Moreover, we show that 
the number of condensate particles at the point of instability decreases with 
increasing temperature, and that for the temperature interval considered, the 
normal part of the gas is stable against density fluctuations at this point. 

PACS numbers: 67.40.-w, 32.80.Pj, 42.50.Vk 
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I. INTRODUCTION 



Since several groups reported evidence for Bose-Einstein condensation (BEC) in atomic 
gas samples of 87 Rb 0, 23 Na |2| and possibly 7 Li || last year, there has been an increased 
interest in this field of physics. After almost two decades, one has finally been able to verify 
the prediction which Bose and in particular Einstein launched for so many years ago. This 
opens a totally new era regarding research on degenerate systems, and it is to be expected 
that there will appear many new and interesting experimental and theoretical results about 
the properties of Bose condensed atomic gases in the near future. 

This is certainly true for the 87 Rb and the 23 Na systems, but maybe less so for the 7 Li 
system, because the experimental results on the latter gas are not yet completely understood. 
In contrast to the former systems where the atoms have a positive s-wave scattering length 
a, 7 Li atoms have a negative s-wave scattering length. This leads to an effectively attractive 
interatomic interaction, which makes the system unstable at large densities. Indeed one of 
us showed that a dilute, homogeneous gas of atoms with a negative s-wave scattering length 
a collapses to a dense (liquid or solid) state before the density is reached at which BEC 
is expected to occur at a given temperature ||. Furthermore, it was shown that the BEC 
transition is actually preceded by a BCS-like transition to a superfluid state and that this 
transition also occurs in the unstable regime of the phase diagram. 

However, as was pointed out by Hulet || and Ruprecht et al. [§] the situation is different 
if the atoms are confined by a magnetic trap. In particular, it was found that for an 
inhomogeneous gas cloud at zero temperature the condensate is stable at the mean-field level 
if the number of particles is sufficiently small, or more precisely if N < NQ >max = 0.573//|a|, 
where I = ^Jh/mu is the typical size of the one-particle ground state in the trap. However, 
quantum fluctuations cause a decay of the condensate on a timescale which is fortunately 
much longer than the timescales at which the experiments were performed. The same holds 
true at the relatively high temperatures of interest experimentally, even though the decay 
is now caused by thermal fluctuations 0. 
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However, in Ref. |7j only the stability of the condensate was discussed. Although this is 
an important first step, it is clearly not sufficient, because we know from the homogeneous 
case that also the non-condensed gas can be unstable against density fluctuations. Therefore 
we will consider here the stability of the complete system including condensate as well as 
above condensate particles. We will try to answer the question when the system as a whole 
becomes unstable and whether it is the condensate part or the non-condensate part which 
causes the instability. 

In this stability analysis we cannot make use of the simple local-density approximation. 
The local-density approximation is only applicable in systems for which the correlation 
length £ (roughly speaking the distance over which the particles influence each other) is 
much smaller than the typical trap size / over which the density changes and the system 
behaves locally homogeneous. However, close to the critical temperature the correlation 
length diverges, and the local-density approximation always breaks down. Nevertheless, 
it is valid at the spinodal point if we satisfy the condition that ksT/hu ^> l/\a\ ^> 1 or 
equivalently N ^> N^ max i.e. the total number N of particles in the gas must be much 
larger then the third power of the maximum number of condensate particles N ^ max . Since 
this amounts to iV ^> 10 9 for the trap parameters of the Rice experiment 0, it is clear that 
on the basis of the local-density approximation we cannot decide if a cloud of 10 4 — 10 5 atoms 
is mechanically stable and therefore (meta)stable. On the other hand, if it were allowed to 
use the local density approximation at the spinodal point, we would immediately conclude 
that a (meta) stable condensate cannot exist. 

To go beyond the local-density approximation, we will present numerical results for the 
free energy of the system at several temperatures. To do so, we first present in Sec. [H] the 
finite temperature theory for the inhomogeneous gas. The equations of motion that describe 
the gas are derived from a variational principle ||. In analogy to the homogeneous case, 
we incorporate the possibility of both a BEC and a BCS transition. Subsequently we give 
an expression for the free energy. Since the experiments with 7 Li are performed at densities 
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and temperatures such that nahq h <C 1 (where n is the density and = y2nh /rnksT 
the thermal wavelength of the atoms) and in particular even at T ~ T c , the theory can be 
simplified by neglecting the possible BCS pairing and using the Hartree-Fock approximation. 
However, in view of the fast experimental developments, it is to be expected that also the 
regime naAZ, 3> 1, which amounts to T <C T c , can be reached in the near future and this is 



why we present here the full theory including also the effects of BCS pairing. In Sec. Ill 



we 



then present our numerical results. In Sec. |III A| we discuss the zero temperature limit of 
the Hartree-Fock approximation and make a comparison with previous calculations. Also a 
thermodynamic criterion for the stability of the gas is given. In Sec. |Tl rB| we proceed to non- 
zero temperatures. We calculate the maximum number of condensate particles as a function 
of temperature and pronounce upon the issue whether the condensate or the non-condensate 



causes the instability of the gas. The paper ends in Sec. IV with some conclusions. 



II. THEORY 

We start this paper with the equilibrium theory for a dilute gas of particles with mass m 
in an external trap potential V ext (r), interacting with each other through an approximately 
local (because \a\ <C I) two-body potential Vo5(f). In the numerical calculations which follow 
subsequently, we will specialize to 7 Li atoms, which have a negative s-wave scattering length 
a and Vq < 0. The interparticle interaction is therefore effectively attractive. 

The grand-canonical Hamiltonian of the system is given by M 

H = J rffj^(f) (- T ^ + V ext (?)-i^j V'(^ + ^W t (r)V' t (r)V'(^(r)|, (i) 

where /i is the chemical potential and ip{r) respectively ip'(r) annihilates and creates a 
particle at position r. As usual, the density of particles in the system n(r) is given by the 
grand canonical average (i/;* (r)i[)(r)) , and the total number of particles is A" = / df n(r), 
which ultimately determines the chemical potential of the gas. 

For particles with a positive s-wave scattering length a, the annihilation operator ijj(r) 
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has a non-vanishing expectation value below the critical temperature. By separating out 
this expectation value in the usual way [[J, i.e. 

V(fO=Vo(r)+V/(r), (2) 

where ip' describes the non-condensate part and if>o(r) = (if)(r)) = yjn (r) is the condensate 
wave function, one can derive the equations of motion for the condensate as well as the non- 
condensate part of the gas. However, in the case of a negative s-wave scattering length, it 
can be shown for the homogeneous case that the free energy as a function of the expectation 
value of the field operator ip has at low temperatures a local maximum for some non-zero 
value of (i/}(r)) and there is only a local minimum for (ip(r)) = 0, which is therefore the 
correct value around which one has to expand ||. As a result we must use a different order 
parameter to describe a phase transition due to quantum degeneracy effects, namely the 
BCS-type order parameter (ij;(f)ip{f}) . In the case of bosons, this is actually known as the 
Evans-Rashid order parameter. 



A. Evans-Rashid transition 

To derive the equations of motion that describe the gas, it is useful to determine the 
(exact) grand-canonical potential 

tt ex (T,n) = -£; B Tln(Trexp (-/?#)), (3) 

from which all thermodynamic quantities that we wish to know can be calculated. It is 
well known that if H t is some trial Hamiltonian, and Q t the corresponding grand-canonical 
potential, we have the variational principle [|Hj 



n ex <n = n t + {H- H t ) t (4) 

where (O)t is the expectation value of the operator O in the grand-canonical ensemble based 
on H t . The trial Hamiltonian that we want to use here is given by 
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Ht 



h 2 v 2 



2m 



+ V ext (r) - n + ft£(r) ip(r) + 



+^A (f)^ t (f)^ t (f) + -A* (r)i[)(r)i;(r)\ . 



(5) 



It is quadratic in the field operators and indeed has non-zero expectation values for both 
(i/j'(r)ip(r))t and {ip(r)ijj(r))t- In this expression, the functions HTi(f), Ao(r) and it's complex 
conjugate A^r) are variational parameters which have to be determined by minimization 
of the grand-canonical potential Q. The trial Hamiltonian has non-diagonal elements pro- 
portional to the BCS order parameter A (r) and the diagonal contribution proportional to 
hH(r) is the self-energy due to the two-body interaction. 
This trial Hamiltonian can be put into the diagonal form 



(6) 



by applying the Bogoliubov transformation 



ijjjr) — (uj(r)bj + v*(r)bj 



(7a) 



(7b) 



The operators and b\ are required to satisfy the usual Bose commutation relations, and 



therefore the functions Uj{r) and Vj(r) are normalized as 



V(#0,V' t (O] = E(' 



5(r — r 1 



(8) 



Using the relations [H t ,bj] = —hujjbj and H t , &j = hujjbj, and substituting Eqs. ( [7a] ) and 
( |Tb"D into the commutators ^ an d H,^ 
solutions to the following eigenvalue equation: 



it is found that Uj(r) and Vj(f) must be 



Hq + H12(r) — \i — hujj 



Anff) 



^ ( u 3 {r) ^ 



V 



/ 



Vj(r) 



(9) 



-A* (f) -H -hE(r) + fi-hojj 

where H = — h 2 \/ 2 /2m + V ext (f). These coupled equations are the Bogoliubov-DeGennes 
equations |TTJ. They can be solved selfconsistently once the functions fr£(r) and A (r) are 



obtained from minimization of Q and expressed in terms of Uj(r) and Vj(r). Furthermore 
the ground-state energy E g in Eq. (^j) is given by 

E 9 = I ^E { -^K'(r)| 2 " ^A oM *(r>,(f) + ±A* (f) Uj (r)v*(r)} . (10) 

We now return to the calculation of the thermodynamic potential Q. Since H t is diagonal 
according to Eq. (P) , it is easily verified that the first term on the right-hand side of Eq. (|J) 
is given by 

n t = E g + fc B rX> (i - e-^), (ii) 

whereas the second term can be rewritten, using Wick's theorem ||, as 

(H - H t ) t = J dr{y (^ t (r)V'(0)t(V' t (OV'(r))t + ~ WW^OMVW^t 

-nS(0(^(r)^(r)) t - -AotrO^rl^r))* - ^A*(r)^(r)^(r))*} • (12) 
Substituting the Bogoliubov transformation from Eq. (|?p, we find that 

<^(rw*o>t = e {(i«,-(#or + ivjimm^j) + i^c^i 2 } , (13) 

and 

(^(r-)^(r)) t = £ Mj (r>;(f)(l + 2iV(^)), (14) 

where the function N(hoUj) = (bjbj) t = [eP hw i — l)^ 1 is the Bose distribution function for 
the Bogoliubov quasiparticles. 

The still unknown functions hH(r), Ajj(r) and Ao(r) have to be chosen such that the 
functional fi[7i£, A , Aq] is minimal, i.e. 

no 

0. (15) 



an 




on 




on 








0A 




dA* 





The last condition is just the complex conjugate of the second, and it suffices to consider 
only one of them. In the Appendix it is shown that Eq. (|15|) requires that 

fc£(f) = 2V Q (V t (r)V(r)> i (16) 
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and 

A„(r) = Voi^^t. (17) 
As is explained in Ref. fll2"l , to incorporate all two-body scattering processes in this many- 



particle system, the factor V in Eq. ( |16D must be replaced by the many-body T-matrix 
T MB , but in the regime of interest where the temperature is large compared to the average 
interaction energy, i.e. naA 2 h -C 1, this can be approximated by the two-body scattering 
matrix T 2B = Aitah 2 /m and we find the usual Hartree-Fock contribution to the self-energy 
ft£(r) = 2n{r)T 2B '. In addition, Eq. (|17D corresponds to the gap equation of BCS theory. 
Since this theory already incorporates all ladder diagrams, the factor Vq should here not 
be replaced by the many-body T-matrix. Collecting together all terms, we find for the 
thermodynamic potential 

llAnff)! 2 " 



Q = E g + k B T^2 H 1 ~ e'^i) - I "df 
j 



n (r)T 



2B 1 l^ol 



2 Vr 



o 



(18) 



From this expression the free energy of the system can be calculated directly using the 
thermodynamic identity F = Q + fiN. 

The equations obtained thus far are only valid when there is no Bose condensate present. 
However, as in the homogeneous case, it is evident that the lowest energy Tlujq will go through 
zero at sufficiently low temperatures and at this point the corresponding one-particle ground 
state becomes macroscopically occupied, i.e. a Bose condensate is formed. Hence, this ground 
state has then to be considered explicitly. 

B. Bose-Einstein condensation 

We now adress the changes in the above equations that are required if a Bose condensate 
is present. First, we consider the limit T — >• 0, for which all particles in the system tend to 
occupy the ground state. In the (trial) grand-canonical ensemble 

Z t = Tr [exp (-f3H t )} = Tr [exp (-(3E g - (3 £ %Ujb%) (19) 



used thus far, one can calculate by standard-statistical physics methods, that for any j the 
expectation value (ftj-fttfyft^ = (Nf) t -(N j ) t = 2(Nj) 2 , where (Nj) t = -(1/ p)d\nZ t /dhuj j = 
N(tiUj). In fact, the factor of 2 in the self-energy Eq. fll6|) originates from this property. 
However, at zero temperature all particles will tend to the ground state and we therefore 
find that the fluctuations in the total number of particles are given by 

(N*) - (N) 2 _ (iV 2 ) - (No) 2 ^ 1 + J_ = 1 + J_ {20) 



(N) 2 (No) 2 (N ) (N) 



and of order O(l) instead of the usual 0{l/y (N)). Hence the fluctuations in the particle 
number are as large as the average itself, which leads to the conclusion that the use of the 
grand-canonical ensemble does not lead to an appropriate description of the Bose condensed 
gas. 

There are basically two ways to restrict the grand-canonical ensemble and circumvent 
this problem. The first one is to introduce a condensate expectation value according to 
) = a//Vq. Clearly, we then have (Nq) = Nq, and at zero temperature we end up with 



hH = 71qT 2B instead of Eq. For non-zero temperatures the effect will be that the 

ground-state wave function satisfies Eq. (Q) with KE = (2n — uq)T 2B , whereas the excited 
state wave functions obey the same equation except that now hH = 2nT 2B . Moreover, the 
free-energy density now contains the term / = \n 2 T 2B instead of n 2 T 2B (c.f. Eq. (pD and 
use / = Q/V + fm). Thus from energy considerations, it is indeed favorable to introduce 
a condensate expectation value if the s-wave scattering length is positive (T 2B > 0), but 
this is not the case for negative a. As mentioned previously, it can be shown that the 
thermodynamic potential is a mexican hat shaped function with extremum at | (bo) \ = 

v^/Vo, which is however inverted with respect to the (bo) -plane when the scattering length 
a changes from positive to negative ||. Therefore, the local minimum at (bo) = v^o that 
is present for a > becomes a local maximum for a < 0, and the local minimum of the 
thermodynamic potential occurs at (bo) = in the latter case. So, for a gas of 7 Li atoms, the 
use of the order parameter (bo) appears not to be the correct way to control the fluctuations 
in the number of condensate particles. 
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The second method to restrict the condensate fluctuations is to introduce a different 
restricted grand-canonical ensemble according to 



Tr [e~^ 
ETr 



-BEa-B V . fiWjbtbj r 



b [ b ,N 



N 

^Q-PhM) (21) 



No 



where 



e -pOt(/j,,No) e -0E g -3hwoNo J~J (_ j n ^ _ Q-P^i^y 

From these expressions it follows that 

iV ) = £ 5 + ^ iV + k B T In (1 - e-^)- 

The total number of particles is therefore given by iV = J2j(Nj)t — 12 j dVL t /dTiujj = N + 
Yj&oN(hujj). Moreover, the largest contribution to the sum over the number of condensate 
particles in Eq. (|2"T|) comes from a minimum in Q t , which implies that Nq = if Tiujq < 
and < N < oo if Tuoq = 0. 

The expectation value (Nq) calculated in this restricted grand-canonical ensemble (by 
construction) equals (N ) 2 , whereas for the other energy levels (j ^ 0) nothing has changed 
compared with the results in the original grand-canonical ensemble. In conclusion we there- 
fore arrive at the same equations for the ground state and the excited state wave functions 
as we had derived by the first method for a gas with positive s-wave scattering length: The 
ground-state wave function (uo(r), Vo(f)) satisfies Eq. (§) with hT, = [2n(r) — n (f)]T 2B 
whereas the excited states have just hH = 2n(r)T 2B . In addition, the condensate den- 
sity obeys no(r) = No(\uo(r)\ 2 + |^o(^*)| 2 ) + I u q(t)| 2 an d the total density is given by 



n r 



n oif) + J2j^o{N(hujj)(\uj(r)\ + \vj(r)\ ) + 1^(^)1 }• The change in expectation 



value (N?) for j = will of course also change the free energy if the system is Bose con- 
densed. Taking this change into account, the grand-canonical potential Q turns out to be 
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given by Eq. fll8|) except that the term j = must be excluded from the summation over 
states and a term | / dr n^{r)T 2B must be added. As a result, the free energy becomes 



F = fi fdf n(r) + E g + A; B T ^ In (1 - e ~^) - f dF 



(n (r) - -n (r))T - - — 

(22) 



Note that the self energy hH for the groundstate contains a term no(r)T 2B . According to 



Ref. [JE^, the two-body scattering matrix in this term should again have been the many- 
body T- matrix. However, at sufficiently high temperatures such that naA 2 h <C 1, or even in 
the opposite regime if the interaction energy between the atoms is smaller than the energy 
splitting of the one-particle states, i.e. nT 2B < %u, the many-body T-matrix T MB can be 
approximated by T 2B . 



C. Mechanical stability 

For the homogeneous gas with effectively attractive interactions, the free-energy density 
satisfies df/dn = fi, and the chemical potential /i as a function of the density n becomes 
multivalued (d 2 f/dn 2 changes sign) for smaller densities than those needed for BEC. This is 
the instability criterium for the homogeneous system. A detailed analysis shows that in this 
case the BEC transition is indeed preceded by a BCS transition, but that both transitions 
occur in the unstable regime of the phase diagram. 

However, in the introduction we already mentioned that in an inhomogeneous system a 
metastable condensate can exist if the number of condensate particles is sufficiently small, 
i.e. Nq < 0.573//|a| 0. Qualitatively this can be understood from the fact that a collapse of 
the condensate requires that other harmonic oscillator states need to be mixed into the wave 
function of the condensate. For this energy is needed (virtually), which can be supplied by 
the interactions provided that the densities in the system are sufficiently high. As a result, 
when the density of the gas becomes so high that the system will collapse, there must be 
some radial unstable mode in the density fluctuations. 
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Systems of compact objects such as white dwarfs and neutron stars |17j can also be 
unstable for collapse under certain conditions. Although a compact object consists of de- 
generate fermions (electrons, protons and neutrons), and furthermore time and length scales 
for stellar systems cannot be compared with those of the Bose condensate, the physics in 
both systems has interesting similarities. 

Indeed, the final state of a star when it has burnt up all it's nuclear fuel is a white dwarf, 
neutron star or black hole, depending on it's mass. Such a compact object is formed because 
in the final stages the star still radiates energy at the expense of gravitational energy, i.e. the 
system contracts. This cannot go on indefinitely, because at a certain point, the electrons 
and protons in the star become degenerate. This causes an extra internal pressure and 
the star will come to equilibrium. For a white dwarf, which has a maximum mass of 1.4 
times the solar mass, this occurs at a radius of about 5000 km. In this stage the gradient 
of the pressure just cancels the gradient in potential energy. However, when the mass of 
the original star is between 1.4 and 3 times the solar mass, the gravitational force will be 
so strong that equilibrium can only be restored when almost all electrons and protons are 
squeezed together to neutrons by inverse /3-decay: In that case the star contracts to an even 
more compact neutron star with a radius of about 10 km. Above these definite maximum 
masses, the white dwarf respectively the neutron star can not support themselves against 
gravitational collapse and this can lead to black hole formation. To study the stability of 
these systems, it is known that it is convenient to parametrize all equilibrium density profiles 
n(r; n c ) by the central density n c of the star and that at the point of instability the mass of 
the object as a function of the central density of matter n c exhibits an extremum. 

In analogy, we thus expect that in the case of a trapped atomic gas, the onset of the 
instability is determined by the condition OF/ dn c = and that there exists a zero mode in 
the density fluctuations at this point. To see this more explicitly, we consider the free energy 
functional F[n], which gives the free energy of the equilibrium density profile n(f; n c ). As a 
result we have 
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F[n(f; n c ) + 5n(f)] = F[n(f; n c )] + j df fi(n c )5n(f) + 

5 2 F 



+ - f dfdf 
2 J 



5n(f)5n(f) + .... (23) 

n(r;n c ) 



8n(f)8n(f 

where fJ,(n c ) = 5F / 5n{f)\ n ^. nc y If the central density is changed slightly, we have 

F[n(f; n c + 8n c )\ = F[n(f; n c ) + 8n c + 0(8n 2 c )} 

= F[n(f;n c )} + J df fi(n c )^^-5n c + 0(5n 2 c ). (24) 

We thus conclude that if dF/dn c = 0, then either /i(n c ) = or / df dn(f;n c )/dn c = 0. 
The latter possibility is in general the physically relevant one because it shows that the two 
density profiles n(r; n c ) and n(r; n c + Sn c ) have the same total number of particles, i.e. 

N{n c + 5n c ) = J df n(f; n c + 5n c ) ~ J dfn(f;n c ) = N(n c ). (25) 

The fact that two density profiles containing the same number of particles have the same 
free energy up to first order indicates that there is a zero mode present. Roughly speaking, 
it does not cost energy to deform the first density profile continuously into the second, which 
indicates the threshold for instability. We therefore anticipate that the onset of instability 
occurs if the free energy has an extremum which conserves particle number (i.e. fi(n c ) ^ 0) 
as a function of the central density of the gas. So, although collapse in compact objects 
respectively in a Bose condensate is caused by a different mechanism, the final criterion in 
both system may be, surprisingly enough, the same. 



D. Hartree-Fock approximation 

We have derived the equations that describe an inhomogeneous gas at non-zero tem- 
peratures. A convenient procedure to solve these equations numerically, would be to start 
with some suitable initial distribution of particles n(f) and an initial BCS order parameter 
Ao(r), and iterate the equations to selfconsistency. It is, however, well known that it is 



rather difficult to ensure the selfconsistency of Ao(r) | IT| , and furthermore that Eq. (|T 
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contains a divergence, because the interparticle potential was approximated by a (^-function 
potential. For a homogeneous gas, this divergence can easily be corrected for, but in this 
inhomogeneous case, it is not a priori clear how one has to deal with it properly, although 
it is not difficult to convince oneself that the divergence can be canceled by calculating the 
molecular states of two atoms in the trap. Fortunately it is not necessary to solve these 
problems here because we are primarily interested in the regime naA 2 h <C 1, where the av- 
erage energy k B T of the particles is much larger than the interaction energy and the effect 
of Ao(r) is very small. Therefore, we neglect in the following the BCS order parameter 
Ao(r), which in turn means that the functions Vj(r) = 0. So, for the uncondensed gas, 
the Bogoliubov-DeGennes equation Eq. (H) then reduces to the Schrodinger equation for a 
particle in an effective potential V e ff(r) = V ext {f) + 2n(r)T 2B 

[ ~ ~2^T + Vext ^ + 2n ^ T2B -Hfa ( f ) = huo ^i f) • (26) 
and the free energy of the system is 

F = n + fiN = k B T ln (! - ex P i-PhWj)) ~ I drn 2 (r)T 2B + fj, f dfn(r), (27) 

j 

where the particle density is given by 

n(r)=J2\H^\ 2 N(hu J3 ), (28) 

3 

and the wave functions <f>j(r) are subject to the condition 

df\(j) j {r)\ 2 = l. (29) 



When the ground state is macroscopically occupied, the noncondensed particles satisfy 
Eq. ( p6|) for j ^ 0, and the condensate wave function satisfies 

("1^" + Vext ^ + [2n ^ ~ n °^ ]T2B ~ Mr) = huJoMf, (30) 

where the condensate density is given by no(f) = N \(po(r)\ 2 and N is the total number 
of particles in the condensate. The non-condensate density is n'(r) = J2j^o Ifijif^Nfiujj), 
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and the total density is n(r) = no(r) + n'(r). The free energy in this case is according to 
Eq. (p2"D given by 



F = k B TY / In (1 - exp (-/3hUj)) - I dr \n 2 (r) - -n 2 (r))T 2B + ll f dr n(r). (31) 
m J 2 J 

Note that for T = 0, i.e. all particles in the groundstate, Eq. fl3"0| ) corresponds to the non- 
linear Schrodinger equation (NLSE), studied for example by Ruprecht et al. || and first 
derived by Goldman et al. ]14| in their pioneering work on spin-polarized atomic hydrogen. 
In addition, note that Bergeman Jl3[ in his analysis uses Eq. (|26| ) with T 2B replaced by 
T 2B jl for both the condensate and the excited state wave functions, which corresponds to 
the Hartree approximation. Although this gives correct results at zero temperature, this is 
no longer true for non-zero temperatures because it does not properly take into account the 
mean-field interactions due to the non-condensate part of the gas. 



III. RESULTS 

At this point we have all tools available to study the stability of atomic 7 Li for temper- 
atures obeying naA 2 h < 1. As mentioned previously, this is done numerically by solving 
Eq. (p6| ) and Eq. (|30|) selfconsistently with the total particle density n{f). The gas is assumed 
to be confined by an isotropic harmonic oscillator potential 




2 2 

mw r , 



where we take for uj the 'average' (uJxUJyUJz) 1 ^ of the (non-isotropic) trap frequencies used 
in the Rice experiment 0. This results in an energy splitting of Tiuj/ks = 7.1nK. Due to 
this simplification the density profile of the gas will depend only on the distance r from the 
center of the trap. The s-wave scattering length of 7 Li is a = — 27.3ao, where ao is the Bohr 
radius flT5| ]. We first consider the case T = and subsequently present results for non-zero 
temperatures. 



15 



A. The case T = 



In this section all particles are considered to be in the condensate, which is the case 
at zero temperature. This has already been subject to extensive research of several other 
papers, see for instance Refs. || and [|nj. For a fixed number of condensate particles Nq, the 
lowest energy eigenvalue and wave function of the Schrodinger type equation (^Up is solved 
by a numerical integration and the density distribution n(r) = no(r), the chemical potential 
fi and the free energy F are determined from this solution. 

In Fig. [I| we plot first of all \i as a function of N Q . If there are only a few particles in 
the condensate, \i is seen to be equal to ^hui, i.e. the groundstate energy of a particle in a 
harmonic oscillator. However, as No increases, the effective potential V e ff(r) = ^muj 2 r 2 + 
n(r)T 2B , grows deeper and deeper in a small range around the center of the trap since T 2B 
is negative, which pulls the particles more and more to the center of the trap. This decreases 
the value of the ground-state energy and consequently also the chemical potential. As can 
be seen from the figure, for Nq > 1241, a solution cannot be found anymore, indicating that 
the condensate becomes unstable. The maximum number of N ^ max = 1241 corresponds well 
with the condition N 0j7nax ~ 0.405 l/\a\ found by Ruprecht et al. for the appropriate trap 
parameters 0. 

The free energy given by Eq. ([TT]) is plotted as a function of N and as a function of 
the central density n c of the gas in Figs. and § respectively. Notice that the derivative 
of the free energy with respect to the number of condensate particles exactly reproduces 
the chemical potential, i.e. \x = dF/dN, showing the consistency of our calculations. From 
Fig. ^| it follows that the free energy as a function of the central density approaches a 
constant (maximum) value. This is also shown in the inset of this figure. So, as anticipated 
in Sec. |II C| , the instability appears at an extremum of the free energy as a function of the 
central density of the gas cloud. If the density increases further, the gas will collapse to a 
dense state. With the theory presented above we clearly cannot describe the gas beyond 
this point, for we would need a theory that can describe the system also at high densities. 
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B. The case T^O 



For non-zero temperatures the particles in the gas occupy the harmonic oscillator states 
in Eq. (|26|) according to the normal Bose distribution at a given chemical potential. However, 
if /i is increased from — oo to some value below ^hu, the number of particles in the ground 
state starts to increase dramatically, and we can only put more particles in the system by 
forcing them into the ground state. At this point, Tluq equals zero and the gas consists 
of a condensate part with density n (r) = A^ o |0o( r )| 2 , as well as a non-condensate part 
with density n'(r) = X^yo N(huj)\(fij(r)\ 2 . So, the calculation of the chemical potential and 
the free energy now consists of two parts which correspond to the use of the unrestricted 
and restricted grand-canonical ensemble respectively. In the first part we increase /i from 
— oo to some maximum value fJ, max , above which there are no longer any solutions. In 
the second part we on the other hand increase the number N of particles in the ground 
state and then determine \x from the ground-state energy of the trap using that hu = 0. 
Note that in this second part the chemical potential decreases again, because the increasing 
density of the condensate lowers the ground-state energy. This behaviour of the chemical 
potential explains why no solutions with /x > Umax could be found in the calculation with 
the unrestricted grand-canonical ensemble. The two parts of the calculation join smoothly 
together within an error of the order of 1% at \i m ax where the condensate fraction is of the 
order of 5%. 

At non-zero temperatures the density profiles of the gas are determined by calculating all 
energy levels and corresponding eigenfunctions up to lOksT. Since the ideal 3D harmonic 
oscillator energy levels are given by e n i = (2n + I + 3/2)huj where n and I are integers, 
this corresponds to taking as many as (1/2 * lOfcsT ' /Tiuj) 2 levels into account. Clearly, this 
number increases rapidly as a function of temperature. When the number of particles in the 
groundstate is small (typically corresponding to fi ~ —Tiuj), Eq. ( |26| ) was used to calculate 
all wave functions and Eq. fl27|) to calculate the free energy. For larger values of fi, the 
groundstate was determined by Eq. ( p0|) and the free energy by Eq. ([H]). To check that 
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our results are consistent, we first compare in Fig. |4] the density profile above the critical 
temperature with the prediction of the local-density approximation, i.e. with 

n(r) = 1/A 3 th g 3/2 (exp [/?(/, - V(r) - 2n(r)T 2B )]), 

The agreement is good for large and negative fi, but for /i ~ — hu, a deviation becomes 
visible around the center of the trap, indicating that the critical temperature is approached. 

In Figs. |] and |] the free energy is plotted as a function of the (total) central density of 
the gas at T = 50nK and T = lOOnK respectively. We again checked that /i = dF/dN as 
required. The inset in both figures shows a magnification of the local minimum in the free 
energy curve where fi goes through zero. This is not a point of instability of the system 
because small variations in the central density do not conserve the total number of particles. 
The point of instability occurs again where the free energy approaches a local maximum. 
Note that the maximum densities in the center of the trap that can be obtained are orders 
of magnitude higher then those for a homogeneous trap. In the homogeneous case, collapse 
occurs already at densities smaller than n BEC = ((3/2)/ A^ h ~ 2.612/A^. For T = 50nK this 
corresponds to ubec — 1-07- 10 11 cm~ 3 , and for T = lOOnK we have ubec — 3.04- 10 11 cm -3 . 
Note furthermore that the interaction term [2n(r) — no(r)]T 2B in the effective potential 
V e ff(r) of the groundstate becomes in the center of the trap as large as n c T 2B ~ —5hu. 

Next we take a closer look at the point of instability of the system when the temperature 
increases. As was pointed out before however, the number of harmonic oscillator states 
which have to be taken into account to calculate the density acurately increases very rapidly 
and thus slows down the calculation considerably. To avoid this, we make for temperatures 
higher than 150nK use of the fact that only the lowest wave functions of the harmonic 
oscillator states are influenced by the interaction term 2n(r)T 2B or [2n(r) — n Q (r)]T 2B for 
the ground state, and the wave functions of the higher states are unaffected, although their 
occupation numbers change, due to the fact that /i equals the ground-state energy if there 
is a Bose condensate present. 

In Fig. [7| we plot at the point of instability the normal density of the gas in the center 
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of the trap n'(0) as a function of temperature (dashed line) and compare this with the 
density n B EC — 2.612/A^ (solid line) required for BEC in the homogeneous case. In the 
inset of the same figure, the number of noncondensed particles as a function of temperature 
is plotted (dashed line), and this is compared with the usual criterion for the onset of 
BEC in a noninteracting gas, i.e. N BEC = ((3)(k B T/huj) 3 ~ 1.202(k B T/huj) 3 (solid line). 
As expected, the consequence of the attractive interatomic interaction is that the non- 
condensed particles are pulled towards the center of the trap. The solid line in Fig. |8| 
shows the maximum number of condensate particles as a function of temperature. Clearly, 
the occupation of the condensate at the point where the gas becomes unstable decreases 
when the temperature increases. In Ref. [|7| it was argued that an increase in temperature 
would not lead to a decrease in the maximum number of condensate particles since the 
non-condensed density is approximately constant over the extend of the condensate wave 
function and therefore only shifts the effective potential V e ff(r) by a constant. If this is true, 
the observed decrease can only be explained by the non-condensed part of the gas becoming 
unstable before the condensate holds the maximum number of particles. This might also 
be physically reasonable because the contribution of the normal part of the gas to the total 
density increases everywhere and especially around the center of the trap. 

In view of this we want to try to answer the question whether it is the condensate or 
(as in the homogeneous case) the non-condensed part of the gas which causes instability. 
Of course, for zero temperature the condensate becomes unstable, and we expect the same 
for such low temperatures that the non-condensed part is only a small fraction of the gas. 
For higher temperatures this might change since the number of non-condensed particles 
increases very rapidly as a function of temperature. To analyse this issue, we calculated in 
the temperature interval < T < 400nK the density profiles of the condensate n Q (r) and the 
normal part of the gas n'(r) at the spinodal point. The condensate fraction N /N decreases 
from 1 for T = to about 0.005 for T = 400nK. Subsequently, we try to add particles to the 
condensate, and try to find a new solution to the non-linear Schrodinger equation Eq. ( |50"D 
for the increased number of condensate particles, while keeping the non- condensate density 
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profile fixed. 

The results are also plotted in Fig. [| The dots in this figure denote the maximum 
number of particles that can be in the condensate given the non-condensate density at the 
temperature of interest. For temperatures T < 50nK, the system becomes already unstable if 
only one particle is added to the condensate, from which we draw the conclusion that at this 
temperature it is still the condensate which renders the instability. For higher temperatures, 
it is possible to add a few particles to the condensate, but this appears te be the result of 
numerical unaccuracies of our calculation. We thus conclude that the condensate is unstable 
at the point of instability of the whole system. 

Because of this result, we suspect that the simple argument that the maximum number of 
condensate particles remains constant because the non-condensed part of the gas is approx- 
imately constant over the extend of the condensate wavefunction, may not be sufficiently 
accurate in this system. The only difference in the condensate wave function at different 
temperatures arises due to the contribution of the term 2n'(r)T 2B to the effective potential. 
In Fig. |] we plotted for T = (and consequently n'(r) = 0) (solid lines), and for T = 300nK 
(dashed lines) the effective potential V e ff(r) = 1 /2mw 2 r 2 + 2n'(r)T 2B and the correspond- 
ing condensate densities n (r). When V e ff(r) for T = 300nK is shifted upward such that 
the zero of both potentials coincide, it is clear that the normal part of the gas effectively 
increases the oscillator strength of the trap potential when the temperature increases. Since 
the maximum condensate size N 0>max oc l/\a\ oc l/y/co, an increase in the effective oscilla- 
tor strength felt by the condensate causes a decrease in the maximum condensate size. A 
measure for the deviation of the effective oscillator strength from the original strength u is 
given by the expression 

2 _ 2 fdr(V ext (r) + 2[n'(r)-n'(0)]T 2B )n (r) 
Weff ^ fdrV ext (r)n (r) 

and for non-zero temperature we thus estimate that the maximum number of condensate 
particles is given by 
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(T) = N 

V U eff(X) 

For T = 300nK this amounts to iVo,ma:r(300) = 1174.5, which corresponds rather accurately 
with the value 1173 from our full calculation. For some other temperatures the maximum 
occupation of the condensate determined in this way is denoted in Fig. [8] by the open circles. 
We can conclude that the growth of the normal part of the gas occurs at the expense of the 
condensate when temperature increases. 



IV. CONCLUSION 

We performed a numerical calculation to study the stability of a Bose condensate in a 
trapped gas of 7 Li atoms at zero and non-zero temperatures. This was done by determining 
all quantum states for particles in a harmonic oscillator trap and interacting via two-body 
scattering. The proposed criterion that the gas becomes mechanically unstable when the free 
energy of the system as a function of the central density of the gas approaches a maximum 
value, is confirmed by the calculations. 

For zero temperature, the maximum number of condensate particles is in agreement with 
previous calculations, and for non-zero temperature this number decreases considerably. 
This is due to the fact that the condensate experiences an effective oscillator strength due 
to the presence of the non-condensed part of the gas. This effective potential increases as 
temperature increases and therefore results in a decrease of the maximum occupation of the 
condensate. For the temperature interval < T < 400 the condensed part of the gas renders 
the instability at the spinodal point, so in contrast to the homogeneous case, the normal 
part of the gas remains stable against density fluctuations. 

Furthermore, from the results in Sec. [jl[B| and the discussion in Sec. fl~B[ it can be 
concluded that at low temperatures it seems necessary to include also many-body effects in 
e.g. the scattering length, since the average interaction nT 2B becomes substantially larger 
than the energy splitting Two. To do so appears to be an important challenge for the future 
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which is not only difficult in practice but even in principle due to the presence of infra-red 
divergences in the theory of the dilute Bose gas Closely related to this issue is the effect 
of the BCS transition on the properties of the gas, which still needs to be incorporated in the 
numerical calculations. Once the experiments enter into this low temperature regime where 
naA^ h 3> 1, it should be interesting to compare the experimental data with the mean- field 
analysis presented here, and to see if possible deviations can be understood by the above 
mentioned corrections. This is of course not only true for 7 Li, but also for any other atomic 
species with a negative scattering length such as 85 Rb and 123 Cs. 
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APPENDIX 

To calculate the minimum of the grand-canonical potential Q of Eq. @, it is convenient 
in the following to introduce a compact notation for the inner product of two states. From 
Eq. (ffl), it is found that the normalization can be rewritten as 



dr (Mr) | 2 - \vj{f)\ 2 ) 



dr 



(ilili), 



-v j [r) 



v 1 ) yv^r) J 



(A.l) 



which thus defines our convention for the inner product. 

We start with the derivative of Q with respect to TiT*. It is given by 
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dtt dE n 



+ E n{k 



dhujA 
) 



+ / dr{2T 2B (^(r)ip{r)) 



d(^(r)j)(r)) t 



2 dKE WW* 5 ))* 



-(^ r (r)^'(r)) i - {ip J {r)ip{r)) t 



0" 



2 o(r) M 



d(ifj(r)tfj(r)) t 



(A.2) 



Equating the whole expression to zero, and grouping together the terms proportional to the 
derivative of an expectation value with respect to KE, the solution is seen to be given by 
Eq. ( |16) and Eq. (0), if we can prove that 



(A.3) 



This is most easily achieved by assuming that all functions Uj(r), Vj(f) and A (f) are 
real. In that case, the first term on the left reduces to 



E 



, _ dTiujj 2 , , 9 



dfv 2 (f) 



(A.4) 



The derivative dTiojj/dTiL can be calculated by perturbing the Hamiltonian according to 



5/7 



^ ^ 



V 







The energy shift in %Uj is then, to first order, given by 



Shojj = (j\5H\j) = / df (u 2 (f) +v](r))6KE, 



and therefore 



dhujj 



df (u 2 (f) + v 2 {fu . 



(A.5) 



Using furthermore that / df (u 2 (f) — v 2 (f)^j = 1, it is found that 



(A.6) 



so Eq. (|A.4 ) can be rewritten as 
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1 / dhujj \ dhujj 
2y \dKE ~ J dKL 



2 V 3 dfiY? 



E / df 



1 



dhujj 



d 2 Tiujj 

— huj 



(A.7) 



The factor between brackets is zero. This can be seen by first applying a Taylor expansion 
to show that the energy change 5hu)j due to the shift KE + SHU can be written as 



Therefore, the last term between the brackets in Eq. ([A.7|) can be identified with twice the 
second-order energy shift. Making then use of the standard expression for this second-order 
energy shift 



i¥=3 



hujj — hui 



(A.8) 



which still holds with our new definition of the inner product, we obtain 

t-^ ( ( dhujj\ 2 , d 2 UujA 



3 \ «7=7 3 1 



S(l-«|(|§) 2 W 



0. 



where the third line follows from the completeness of the eigenstates of the Bogoliubov- 
DeGennes equation and from the fact that if we write 2hujj = (hujj — Tiuji) + (tiUj + huii) only 
the antisymmetric part contributes. Finally, we also need that 

d o\ 

-1 



dH 



\ 



J 



Collecting all terms together we thus indeed find that (see Eq. (|A.3| )) 



3 



J dfJ2 {v](r) + N{%Uj){u)(T) + v 2 {f))} - J df (^(f)^(f)) 



0. 



(A.9) 
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At the same time the derivative of with respect to Ao must be zero. That this is also 
the case can be shown by a similar calculation. Assuming A , Uj and Vj to be real, the 
equation dQ/dA = reduces to 



dE g 

OAr 



if the functions frE and A are given by Eqs. 
on the left-hand side is equal to 



(A.10) 
and (|17D . Using Eq. ( |A.6| ), the first term 



dEg 

dAr 



dAr 



df v 2 Af) 



1 v-^ f dhuj ( dUujj . N 



OA, \ dKL 



(All) 



As in the previous calculation, the first and second derivatives of Tiujj with respect to Ao 
and frE can be calculated by perturbing the Hamiltonian according to H = H + 5H, with 



5H 



( 5KL 5A ^ 



\ 



-SA Q -5hE 

The energy levels then shift to second order according to 



dA ( 



+ ; 



( S^E 2 + 2 f^i SfiESAr 



2 OLA Q 



dAl 



(A.12) 



2 \dhZ 2 ~" ' dhE8A 

Comparing this again with the perturbative expression for the energy shift Eq. (|A.8| ), we 
find that 



2 / df Uj(r)vj(r). 



dtujj _ (j\6H\j) _ dH 
dA 5A Ul dA U) 

Moreover, the second order derivative in Eq. (|A.11|) can be written as 



(A.13) 



d 2 hujj 
dA n dhE 



E 



{j\mmWz\j)+{j\Wz\mm-„\j) 



8H 
I dhY, 



8H 
dhY: I 



8H 
I 8A 



Therefore Eq. ( |A.ll ) becomes 



%ujj — huji 



(A. 14) 
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(A.15) 



since 



«9F 



f 1^ 



V 



-1 



and 



dH 

Ms 



-1 



Combining all results, Eq. ( JA.10 ) becomes 



£ / ^(f)^(f)(l + 2A(^)) - i(^) t - ^ty f )t = 



(A.16) 



which completes the proof. 

To gain even more confidence in our expression for Q, it is useful to see whether 

_ dfl 
dfi 

is satisfied. To prove this, we first notice that n 2 {r) = (fc£(r)/2T 2S ) 2 and that according to 
the matrix form of the Hamiltonian (cf. Eq. @) 



N 



J dr n{r) = J df {ip\r)ip{f)) t 



(A.17) 



dhujj 



dhujn 



d^ dhZ ' 

Therefore, we can make use of our previous results in Eqs. ( |A.5|) and (|A.7Q to obtain 

dQ dE„ AT/ . .dhujj 

7T = IT- + N ( huJ j)^ 



(A.18) 



dfx 
jdE : 



9 + Nihujj 



dhujj 



\dhE ' ~' v " " dUE t 
E / & {(\uj(r)\ 2 + IvjWNQiu,,) + |^-(r)| 2 } 



— / drn(r), 



as desired. 
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FIGURES 



FIG. 1. Chemical potential as a function of total number of condensate particles. 

FIG. 2. Free energy as a function of total number of condensate particles. 

FIG. 3. Free energy as a function of the central density. The inset shows that the derivative of 
the free energy with respect to the central density approaches zero at the point of instability and 
that at this point the number of particles as a function of the central density exhibits a maximum. 

FIG. 4. Comparison between exact density and 1/A| h g 3 / 2 (C) above the critical temperature 
for (1) fi = —huj, (2) fi = —lOhto and (3) fi = —2b%uj. 

FIG. 5. Free energy as a function of the central density for T=50nK. 

FIG. 6. Free energy as a function of the central density for T=100nK. 

FIG. 7. Normal density n'(0) in the center of the trap at the point of instability of the sys- 
tem (dashed line). The solid line shows ubec = 2.612/Aj^. The inset shows the number of 
non-condensed particles as a function of temperature for the harmonic oscillator with (dashed) 
and without interaction (solid), which is given by N' = 1.202(/cbT/?uj) 3 

FIG. 8. Maximum number of condensate particles (solid line). The dots denote the maximum 
condensate number at fixed normal part, and the open circles give the maximum occupation of the 
condensate that can be calculated from the effective oscillator strength due to the presence of the 
non-condensed part of the gas. 

FIG. 9. Condensate density (right scale) and effective potential (left scale) for T = (solid 
lines) and T = 300nK (dashed lines). 
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